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ABSTRACT 


A  numerical  study  was  performed  to  investigate  the  dynamic 
crack  propagation  in  fibrous  composite  plates  utilizing  the 
finite  element  method.  A  rectangular  plate  of  uniform 
thickness,  which  had  a  propagating  central  crack,  was  used  for 
the  study.  The  plate  was  a  unidirectional  composite  panel  and 
the  load  was  applied  in  the  longitudinal  direction  of  the 
composite  plate.  Fracture  energies  were  calculated  for  given 
speeds  of  cracks.  The  objective  of  the  study  was  to  examine 
the  effect  of  different  composite  material  properties,  crack 
speeds,  and  densities  on  the  fracture  energy.  The 
discontinuous  node- release  technique  was  used  to  model  the 
crack  propagation. 

The  numerical  study  showed  that  the  fracture  energy  was 
higher  for  a  lower  elastic  modulus  ratio,  if  all  other 
conditions  were  held  the  same.  Furthermore,  a  lower  crack 
propagation  velocity  or  a  lower  material  density  resulted  in 
a  higher  fracture  energy,  respectively,  provided  the  rest  of 
the  parameters  held  constant. 
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I. 


INTRODUCTION 


Dynamic  fracture  deals  with  problems  having  a  crack 
within  a  body  where  the  inertial  forces  play  an  important 
role.  The  early  investigation  of  the  dynamic  fracture  analysis 
goes  back  to  the  case  of  unstable  crack  propagations  in  some 
ships  during  World  War  II.  The  dynamic  fracture  problems 
involve  crack  initiation  and  crack  propagation.  There  are  two 
principal  types  of  problems.  One  is  applying  a  dynamic  load  to 
a  body  containing  stationary  cracks,  and  the  other  is 
concerned  about  a  body  with  a  moving  crack.  The  second  type  of 
problems  can  be  classified  into  two  different  modes; 
propagation  and  generation  modes. 

Stationary  crack  problems  involve  cracks  within  structures 
that  are  subjected  to  dynamic  loading,  such  as  an  impact  load. 
For  example,  Chen  performed  a  numerical  analysis  to  find  the 
dynamic  stress  intensity  factor  for  a  rectangular  plate  having 
a  central  crack.  He  used  the  finite  difference  technique  and 
discussed  the  dynamic  behavior  of  the  stress  intensity  factor. 
He  argued  that  the  oscillations  on  the  stress  intensity  versus 
time  curve  were  mainly  due  to  scattering  phenomena  from  the 
crack  tip  and  the  boundary  surfaces.  He  concluded  that  the 
overshooting  of  the  dynamic  stress  intensity  factor  is 
attributed  mainly  to  the  geometry  and  loading. [Ref.  1] 
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The  generation  mode  involves  moving  cracks  expanding  at 
prescribed  rates.  The  main  objective  of  the  study  is  to 
compute  the  stress  intensity  factor  or  the  fracture  energy  and 
to  determine  the  fracture  toughness.  The  propagation  mode 
deals  with  the  cases  where  fracture  toughness  is  known.  The 
objective  is  to  predict  the  crack  propagation  and  its  arrest 
for  a  given  load.  One  of  useful  models  to  investigate  the 
crack  propagation  and  arrest  was  the  double  cantilever  beam. 
Kanninen  and  his  coworkers  [Ref.  2]  proposed  the  model  and 
studied  it  extensively.  Some  of  studies  in  the  dynamic 
fracture  were  listed  in  the  references  [Ref.  3-7]. 

Although  there  are  many  literatures  available  for  the 
dynamic  fracture  of  isotropic  materials,  the  studies  of 
composite  materials  are  not  so  common.  Some  of  the  fracture 
studies  of  composites  were  shown  in  references  [Ref. 8 -10]. 

This  study  considered  a  rectangular  plate  made  of 
unidirectional  composite  and  containing  a  central,  moving 
crack  subjected  to  a  uniform  tensile  load.  The  finite  element 
method  utilizing  bilinear  rectangular  elements  was  used  for 
this  study.  A  parametric  study  to  compare  fracture  energies  in 
composite  plates  for  various  conditions  was  conducted. 
Different  elastic  modulus  ratios,  crack  velocities,  and 
material  densities  were  considered. 
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II .  MATHEMATICAL  DERIVATION 


A.  GOVERNING  EQUATIONS 

The  governing  equations  for  the  problem  of  elasticity  are 
equations  of  equilibrium,  constitutive  equations,  and 
kinematic  equations.  For  the  two-dimensional  problem,  they  are 
given  below.  Equations  of  equilibrium  are 


d°X^  dxxy_  &U 

dx  dy  P  0t2 

dx*y  +aqy_D 
dx  dy  v  dt2 


=0 

=0 


(1) 


where  ax,  ay,  and  are  the  two  normal  stresses  in  the  x  and 
y  direction  and  the  shear  stress,  respectively,  p  is  the 
material  density,  and  u  and  v  are  the  displacements  in  the  x 
and  y  direction,  respectively.  In  addition,  x  and  y  are 
rectangular  coordinate  axes  and  t  denotes  time. 

Constitutive  equations  are  also  known  as  stress -strain 
relations.  The  relations  are  given  in  a  matrix  form  like 


{o}=  [D]  (e)  (2) 

where  {a}  and  {e}  are,  respectively,  stress  and  strain 
vectors,  that  is 

{oMox  ay  {3) 

{e)=kx  ey  Y*y}T 

t 
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and  [D]  is  the  material  property  matrix  of  size  3x3 .  For  an 
isotropic  material,  the  [D]  matrices  for  plane  stress  and 
plane  strain  conditions,  respectively,  are 

l  u  0 
u  1  0 

0  0  — 

2 

and 


t  D] 


011 

012 

013 

021 

022 

023 

031 

032 

033 

(6) 


The  details  of  Eq. (6)  are  shown  in  ref erences . [Ref .  11] 

Kinematic  equations  are  also  known  as  strain- displacement 
relations.  The  infinitesimal  strains  are  expressed  in  terms 
of  displacements  as  shown  below: 
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(7) 


du 
dx 
dv 
dy 

du  +  dv 
dy  dx 

B.  BOUNDARY  CONDITIONS 

On  the  essential  boundaries  displacements  u  and/or  v  are 
prescribed,  and  on  the  nonessential  or  natural  boundaries 
tractions  are  prescribed.  Usually  on  a  part  of  the  ooundary 
tractions  are  prescribed,  while  on  the  remainder  of  the 
boundary  displacements  are  known.  The  tractions  are  given  by 


^ x  a  x^x  +  1c  xy^v 
^y~^xy^x+<Jy^y 


(8) 


where  S>x  and  $y  are  the  tractions  in  x  and  y  directions, 
respectively,  and  nx  and  iiy  are,  respectively,  the  components 
of  outward  unit  normal  vector  to  the  boundary. 

Applying  the  method  of  weighted  residual  technique  to 
the  equations  of  equilibrium  yields 


+  d°X^d'Xy 
dx  dy 

^  ^*y^d°y 
dx  dy 


)  w1d£)=0 
)  c o2dQ=0 


(9) 


where  and  u2  are  test  functions  and  Q  is  the  domain  of  a 
given  problem.  Applying  the  weak  formulation  results  in 
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Jq 


g2u  da1 

P - ~al+Ox-~+X 

dt2  dx 


Cp"V 

p  — — —  0),  +T 

dt 


da. 


0G>11. 

^  ?y  IdQ  +  Jrtxnx°>1+Tjfrnyt,>4dr  = 
^2[ 


+a. 


2  “2  Vxy  0x  ~y  0y 


0  (10) 


and  T  is  the  boundary  of  domain  Q.  Rewriting  Eq. (10)  gives 


Jq(p 


a1  0 

0  G) 


2j 


d2u 


dt2 

d?v 


dt< 


>+ 


da  x 

0X 


0G)1 


0y 

00), 


da. 


dy  dx 


(<r. 


'xy. 


y)  dQ  =  Jrfr“4dT  (ID 


,rl^yt02 


where 


Jq 


d2u 

JoP 

g>x  0 

0  g>2 

dt2 

&v 

da 

dt2 

'  dax 
dx 

0 

da1 

dy 

\°A 

0 

da2 

dy 

da2 

dx 

yy) 

w 

dft 


(12) 


(13) 


and 


(14) 


are  the  mass  matrix,  stiffness  matrix,  and  load  vector, 
respectively. 
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A  finite  element  discretization  is  applied  to  the  domain 
and  the  following  matrices  are  computed  at  the  element  domain. 


C.  HASS  MATRIX 

The  mass  matrix  is  given  by 


Jo*  P 


&U 

0 

dt2 

0  W2. 

dt2. 

dQ 


(15) 


In  this  study  both  bilinear  rectangular  and  linear  triangular 
elements  will  be  used  to  produce  the  finite  element  mesh. 
First  the  finite  element  derivation  will  be  shown  in  detail 
using  the  bilinear  rectangular  element.  The  rectangular 
element  is  shown  in  Figure  1. 


H».c) 

y,v 

IM 

3 

2 

X.U 

4 

1 

(b.cj 

Figure  1  Bilinear  Rectangular  Element 
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The  shape  functions  are 


#i= 
*2  = 
*3  = 


1 

4bc 

1 

4  be 
1 

4  be 
1 

4  he 


(Jb-x)  (c-y) 
(b+x)  (c-y) 
(Jb+x)  (c+y) 
(b-x)  (c+y) 


(16) 


where  b  and  c  are  greater  than  zero.  The  displacements  in  the 
x  and  y  directions,  respectively  are  interpolated  as  shown 
below 


u =HX  ux  +H2  u2  +H2  u3  +Ha  u4 
v=Hx  V,  +H2  v2  +tf3  v3  +Ht  v4 


(17) 


where 


ui 

u2 

U3 

V3 

V-4 


(18) 


is  the  displacement  vector  at  the  nodes  of  the  element. 

Because  the  shape  functions  are  independent  of  time  the 
acceleration  terms  become  as  shown  below: 
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BPu 

dt2 

&v 

dt2 


> 


riUi+H2U2+H^+HiuA 


0  H2  0  H3  0 
0  Hx  0  H2  0  H2 


ha  o 

0  tf4 


(19) 


Substituting  Eq. (19)  back  into  Eq. (15)  and  applying 
Galerkin's  method  results  in 


Ja«P 


*1 

0 

0 

*2 

0 

0 

*2 

0 

0 

*3 

#4 

0 

0 

*4. 

o  h2  0  h3  0 

Hx  0  H2  0  tf3 


*4 

0 


(20) 


Finally  Eq. (20)  becomes 
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Jq»P 


h\ 


Hi  0  HxH2 
0  Hi 
HxH2  0 
0  HXH^  0 
HxH3  0 

0  hxh3  0 

HXH<  0  H2H, 
0  HxHt  0 


H2H3 


HxH3 


Hi 


H2H3 


h2h4 


H, 

0 


h3h. 


0 


0  HxH2  0  hxh3 

0  H2H3  0 

0  H2H3  0 
2 


h2h 4 


0  Hi 

H3Ht  0 


H\H^ 


h2h. 


0  tf3tf4  0 


Hi  0  H3H4 


H: 


dQ 


(u. 


u. 


u. 


u„ 


(21) 


Carrying  out  the  integrations  gives  symmetric  consistent 
mass  matrix  as  shown  below: 


[Mc]  =^bct) 


0 

4 


0  2 
1  0 


0 

2 

0 

4 


(22) 


where  t  is  the  thickness  of  the  element.  Equation  (22)  is 
called  the  consistent  mass  matrix. 

Another  choice  for  the  mass  matrix  is  the  lumped  mass 
matrix.  The  lumped  mass  matrix  is  a  diagonal  matrix  so  that  it 
reduces  computation  time  and  space  needed  for  storage.  The 
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lumped  mass  matrix  is 


10  0 
1  0 
1 


[Afj]  =p  bct\ 


0  0  0 
0  0  0 
0  0  0 
10  0 
1  0 
1 


0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
1  0 
1 


(23) 


D .  STIFFNESS  MATRIX 

The  stiffness  matrix  is  given  by 


3<o  1 
dx 


0 


0 

3g>2 

dy 


3(0^ 

dy 

3g>2 

V\ 

dx 

?xy} 

(24) 


Substituting  constitutive  equations,  Eq. (2) ,  into  Eq. (24) 
yields 


dct]^ 

dx 

0 

3(0,' 

dy 

[D]< 

[M 

0 

3o>2 

3g>2 

ey> 

D f*y. 

dy 

dx : 

dQ 


(25) 
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Using  kinematic  equations,  Eq.  (25)  becomes 


Jo* 


'  3w1 

0 

dti>x 

du 

dx 

dx 

dy 

[D]- 

dv 

do. 

3<o, 

dy 

0 

2 

dy 

2 

dx . 

dv  ^  du 
dx  dy 

dd 


(26) 


The  kinematic  equation  can  be  written  in  terms  of  the  nodal 
displacements  as  below: 


du 

dx 

0 

d 

dv 

dx 

dy 

dy 

>  — 

0 

a 

d 

dv+  du 

dy 

dx\ 

dx  dy 

_a_ 

0 

ai 

dx 

dy 

n 

_a_ 

a 

dy 

dx 

$ 


Hx  0  H2  0  H3  0  Ht  0 
0  Hx  0  H2  0  H3  0  Ht 


u. 


u. 


u. 


u, 


'  dHx 

o 

dH2 

0 

dH3 

0 

dH4 

0 

ui 

V, 

dx 

dx 

dx 

dx 

1 

u2 

0 

an, 

0 

an2 

0 

a#3 

0 

a//4 

* 

ay 

ay 

ay 

ay 

U3 

dHx 

dHx 

dH2 

dH2 

dHz 

dH2 

a#4 

an4 

^3 

]  I 

dy 

dx 

dy 

dx 

dy 

dx 

ay 

ax. 

U4 

V4 

(27) 
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where  the  last  matrix  is  called  [B]  matrix  and  the  last  vector 
called  the  displacement  vector  {u} .  Substituting  back  Eq. (27) 
into  Eq. (25)  yields 


3w1 


0 


0 


002 

dy 


~dy 

dx 


[Z3]  [£]  dQU  . 


(28) 


Using  Galerkin's  method  the  first  matrix  in  Eq. (28)  becomes 


dx 

0 


0 

dy 


'  0tf x 

0 

dH 1‘ 

dx 

dy 

0 

dHx 

dHx 

~dy 

dx 

dH, 

0 

dH2 

dx 

dy 

r\ 

dH2 

dH2 

dy 

U 

dy 

dx 

00)  2 

dH3 

n 

dH3 

0X  . 

dx 

U 

dy 

0 

dH3 

dH3 

~dy 

dx 

dHt 

0 

0H4 

0X 

dy 

0 

0H4 

8H4 

dy 

dx 

(29) 


which  is  equal  to  the  transpose  of  B  matrix. 
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Finally,  the  stiffness  matrix  is 


m  =/Q.  [B]T[D]  [5]dQ.  (30) 

Computing  the  [K]  matrix  using  bilinear  rectangular  element 
gives 

/fcJ^[S]£[D]  [fl]  tdxdy  (31) 

where  t  is  the  thickness,  and 


( c-y)  0  ( c-y)  0  ( c+y)  0  -  ( c+y)  0 

[B]  =  — 0  -(b-x)  0  -(b+x)  0  (b+x)  0  (b-x)  qj) 

C[-  (b-x)  -(c-y)  -(b+x)  (c-y)  (b+x)  (c+y)  (b-x)  (c+y)j 


For  a  general  material  property  matrix  [D]  of  Eq. (6),  the 
element  stiffness  matrix  becomes 


*11 

*12 

*13 

*14 

*15 

*16 

*17 

*18 

*21 

*22 

*23 

*24 

*25 

*26 

*27 

*28 

*31 

*32 

*33 

*34 

*35 

*36 

*37 

*38 

*41 

*42 

*43 

*44 

*45 

*46 

*47 

*48 

*51 

*52 

*53 

*54 

*55 

*56 

*57 

*58 

*61 

*62 

*63 

*64 

*65 

^6 

*^8 

*71 

*72 

*73 

*74 

*75 

*76 

*77 

*78 

*81 

*82 

*83 

*84 

*85 

*86 

*87 

*88 

(33) 


The  detail  of  KAj  is  shown  in  the  Appendix. 


If  a  material  is  isotropic  and  the  plane  strain  condition 
applies,  the  stiffness  matrix  is 

[* a  =  [jg  +  [*,]  .  04) 

[KJ  is  given  by 


~b3c  4 ab2c2  -|-b3c 

-4 ab2c2 

-j^b3c 

-4 ab2c2 

~16b3c 

4 ab2c2 

4^-bc3  lab2c2 

"^6bc3 

-4ab2c2 

—be3 

3 

-lab2c2 

|  be3 

^-b3c 

3 

-4ab2c2 

~\Sb3c 

-4 ab2c2 

- j-b3c 

4ab2c2 

^-bc3 

Xab2c2 

-|bc3 

4ab2c2 

- —be 3 

^b3c 

3 

4  ab2c2 

-|b3c 

-4 ab2c: 

^f-bc3 

4  ab2c2 

~16bc: 

3 

^b3c 

-4 ab2c: 

Jfbc3 

(35) 


where 

A=  E(l-\i)  t 

{Abe) 2  (1+u)  (l-2u) 

a  =  — — — — — 

(1-u) 


(36) 


and  [KB]  is  given  by 
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-y--b3c  4 b2c2 

lb3c 

3 

- 4b2c 2 

^fb3c 

-4 b2c2 

~fb3c 

4b2c2 

If  be3 

4  b2c2 

'fbc3 

-4b2c2 

be 3 

-4b2c2 

|  be3 

±fb3c 

-4  b2c2 

-\eb'c 

-4£>2c2 

zft3c 

4b2c2 

±fbc3 

4b2c2 

be3 

4b2c2 

- fbc 3 

Ifb’c 

4b2c2 

~b3c 

-4b2c2 

Ifbc3 

4b2c 2 

-16 

3  bC 

±fb3c 

-4  b2c2 

be3 

where 


B= 


Et 

2  (4 be) 2  (1+u) 


E.  LOAD  VECTOR 

The  load  vector  is  given  by 


(38) 


Linear  shape  functions  can  be  used  to  find  nodal  load  values 
due  to  the  traction.  This  results  in  lumping  tractions  into 
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two  nodal  points  associated  with  the  boundary  of  the  element. 
If  there  is  no  traction,  the  load  vector  is  zero  at  the 
associated  nodal  points.  Linear  shape  functions  are 


s,*  t-s 


s-s. 


hi=Si*l  Si 


(40) 


where  s  is  the  coordinate  axis  along  the  boundary,  and  is 
the  length  of  the  boundary  of  the  element.  Combining  all  of 
the  above  procedure  gives  an  equation  of  matrix, 


[Af]  (u}+  [K]  {u}=(i?}  (41) 

which  describes  the  equation  of  motion  with  no  damping. 
Damping  was  neglected  in  Eq. (41) . 

F.  LINEAR  TRIANGULAR  ELEMENT 

Since  this  study  deals  also  with  the  linear  triangular 
element,  the  derivation  of  the  matrices  will  be  shown  here, 
but  not  in  detail  since  the  detail  was  shown  before  in  the 
bilinear  rectangular  element  case.  The  linear  triangular 
element,  as  shown  in  Figure  2,  has  three  nodal  points  which 
have  coordinate  values  (xi»yi)  »  (x2,y2)  »  and  (x3,y3) 
respectively.  Linear  triangular  element  has  two  degrees  of 
freedom  per  node. 
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Figure  2  Linear  triangular  element 


Note  that  three  nodal  points  are  numbered  in  the  counter¬ 
clockwise  direction  from  an  arbitrarily  selected  corner. 

The  shape  functions  are 


*i=  za  [Ai+six+ci  n 

H2  =  ±[A2+B2x+C2y]  (42) 

H2  =  ~h  ^3  +S3x+C3y] 


where  A  is  the  area  of  the  triangle,  in  other  words 


and 


det 


1  Yi 

i  *3  y  3 


(43) 
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Ai=x2y3-x3y2  si=y2-y3 

A2  =x3y,  -x3y -j  B2  =y3-y1 

A3=Xiy2-x2yi  B2^y1-y2 


Cr=x3~x2 

C2=Xl~X2 

C2=X2~X1 


It  can  be  noted  that 


A=(A1+A2+Ai)  /2  . 


(44) 


(45) 


The  displacements  in  the  x  and  y  directions  are  given  by 


u=H1u1+H2  u2  +H2  u3 
V~HX  v1  +H2  v2  +H3  v3  . 


(46) 


Carrying  out  similar  calculations  and  procedure  gives  the 
symmetric  consistent,  and  lumped  mass  matrices,  respectively 


[Afc]  ^^-At 
c  4 


0 


0 

1 

0 

2 


(47) 


[Mjl  =^At 


0  0  0  0  0 
1  0  0  0  0 
10  0  0 
10  0 
1  0 
1 


(48) 


The  stiffness  matrix  is  obtained  from  a  combination  of  two 
components 
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! i  BXCXV  BXB2  BxC2 U  BxB2  BxC2v 

C\  C1B2v  C1C2  C1B3v  C1Ci 

B2  B2C2x>  B2B3  B2C3u 

C2  C2fl3u  C2C3 

Bl  B3C3U 

C32 


and 

c2  qs,  c,c2  cxb2  qc3  qB3 

•®1  -®1^2  -®l-®2  -®1^3  -®l-®3 

c22  c,b2  c2c3  c2b3 

b|  b2c3  b2b3 

C32  c3b3 
B2 

where  V  is  the  volume  of  the  element,  that  is, 

V=At  (51) 

and  t  is  the  thickness  of  the  element. 

6.  TRANSIENT  ANALYSIS 

In  this  study  the  central  difference  method  is  used  as  a 
time  integration  scheme.  The  equation  of  motion  at  time  t  is 
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(52) 


The  acceleration  at  time  t  is  given  by 

{u}c=_L.  [{u)t-*t_2{u}t+{u}t*At]  ,  (53) 

At2 

The  velocity  at  time  t  is 

{u}£=_JL  [-{uK^Mu}^]  .  (54) 

2At 

The  error  associated  with  the  acceleration  and  velocity 
calculations  are  of  order  (At2) .  Substituting  the  relations 
for  acceleration  and  velocity  into  the  equation  of  motion, 


( — —  [M]  +  [C]  )  {u)t+4t= 

Ufc)2  2At 

(i?K-  (  [K\  1_  [M] )  {u}c-  ( -L-  [M\  --JL  [C]  )  {u}t_4t  . 

(At)2  (At)2  2a  t 


(55) 


The  method  requires  the  solution  at  time  -At.  Using  the 
equations  of  acceleration  and  velocity,  letting  time  to  be 
equal  to  zero  and  eliminating  {u}At  yield 


-a  t(u}° + {u)o 
2 


(56) 


Note  that  {u}°  can  be  found  using  the  equation  of  motion,  Eq. 
(52)  ,  at  time  zero  since  {u}°  and  {u}°  are  known. 

The  central  difference  scheme  is  conditionally  stable  and 
requires  that  the  time  step  At  should  be  smaller  than  a 
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critical  time  step  Atcr.  The  critical  time  step  is  given  by 

At,  (57) 

where  Tmin  is  the  smallest  period. 

Therefore  the  method  can  be  summarized  as; 

1.  Compute  mass,  damping,  and  stiffness  matrices. 

2.  Compute  {u}°  from 

[Af]  [<T]  {u}°i-  [J^j  (u}°  (58) 

where  {u}°,and  {u}°  are  initial  conditions. 

3 .  Compute 

{u}-At={u)°-A  .  (59) 

2 

4.  Compute  effective  mass  matrix  from, 

[M]  =— -i—  [M]  [C]  .  (60) 

(At)  *At 

5.  Compute  the  effective  load  vector  at  time  t, 

( [ K\  [M\  )Ufc-(— [M]  --JL  [C]  )(u}c-At.  (61) 

(At)2  (At)2  2At 

6. Solve  for  displacements  at  time  t+At  from, 

lfi\  .  (62) 

7. Compute  accelerations  and  velocities  at  time  t, 
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(63) 


(u}c= — i —  [{u}c'aC-2{u}c+{u}c*aC] 

At2 

{u}t= — L_  [-{u}c'AC+{ii)t>4C] 

2  At 

Note  that  steps  5-7  have  to  be  performed  repeatedly  for  each 
time  increment. 

H.  DYNAMIC  FRACTURE  ANALYSIS 

There  are  many  parameters  to  describe  the  crack  tip 
behavior  for  determination  of  the  crack  propagation  and 
arrest.  One  of  the  parameters  is  the  stress  intensity  factor, 
K,  at  a  crack  tip  which  is  defined  as 


where 


K=limz^2nran{r,  0,  t) 


(64) 


oa=on(r,B,t)  (65) 

is  the  stress  component,  and  is  evaluated  near  the  crack  tip. 
In  this  case,  crack  propagation  is  governed  by  stress 
intensity  factor  K,  and  the  material  property  Kc  ,  called  the 
fracture  toughness.  A  crack  becomes  unstable  if  the  stress 
intensity  factor  reaches  the  fracture  toughness  of  the 
material . 
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The  dynamic  analysis  requires  K  be  considered  in  two 


parts.  For  crack  initiation, 

K(a,o ,  t)  =KD(a)  (66) 

where  a  is  the  crack  length,  a  is  the  applied  stress,  t  is 
time,  and  a  represents  the  loading  rate.  For  a  propagating 
crack, 


K(a,o,t)=Kd(&)  (67) 

where  a  indicates  the  crack  speed. [Ref.  7] 

An  alternative  criterion  for  crack  motion  is, 

G(a,  or,  t)  =R(a)  (68) 

where  G  is  the  dynamic  energy  release  rate  and  R  is  the  energy 
dissipation  rate  required  for  crack  growth  which  is  the 
critical  value  of  G. 

The  driving  force  or  the  dynamic  energy  release  rate  for 
crack  propagation  has  contributions  from  strain  energy, 
kinetic  energy,  and  work  done  on  the  body.  The  driving  force 
is  the  net  change  in  these  components  per  unit  area  of  crack 
extension.  The  equation  for  G  is. 
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(69) 


G=1  {  dW-  dU -  dT ) 

baa  da  da 
1  (  dW_  dU  _  dT > 
b&  dt  dt  dt 

where  U  is  the  strain  energy,  T  is  the  kinetic  energy,  W  is 
the  work  done  by  the  external  loading,  a  is  the  crack  length, 
and  b  is  the  plate  thickness  at  the  crack  tip. 

The  dynamic  energy  release  rate  can  be  directly  connected 
to  the  dynamic  stress  intensity  factor.  For  plane  strain 
conditions,  this  relation  is 


G=±^-A(*)K2  (70) 

E 

where  E  and  v,  as  usual,  are  the  elastic  modulus  and  Poisson's 
ratio,  respectively,  and  A  is  a  geometry  independent  function 
of  the  crack  speed  given  by 


A=- 


C2  Cf 


(1-U)  [4  (1--^)  1  (1--^)  2-{2-^)z] 

cl  cl  cl 


(71) 


where  Cx  and  C2  ,  respectively,  are  longitudinal  and  shear 
(transverse)  wave  speeds  in  the  material.  Wave  speeds  are 
given  by, 


(^+4Gs/3) 


1/2 


P 


(72) 
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and 


C2=(Gs/p)1/2 


(73) 


where  Kb  and  Gs  are  bulk  and  shear  modulus,  respectively. 
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III. 


RESULTS  AND  DISCUSSION 


A.  VERIFICATION  PROBLEMS 

In  order  to  verify  the  developed  computer  code,  two 
example  problems  with  known  solutions  were  analyzed  for 
isotropic  plates.  One  was  a  stationary  crack  problem  subjected 
to  a  rapidly  varying  dynamic  load,  and  other  was  a  propagating 
crack  problem  with  a  constant  speed  and  subjected  to  a  uniform 
load. 

The  first  verification  problem  was  a  stationary  crack 
problem  known  as  Chen's  problem  [Ref .1] .  A  rectangular  bar 
with  a  centrally  located  crack,  shown  in  Figure  3,  was 
considered.  A  uniform  tension  of  magnitude  0.004  Mbar  with 
Heaviside- function  time  dependence  was  applied.  The  material 
properties  were  given  as  below;  shear  modulus  of  0.76923  Mbar, 
bulk  modulus  of  1.66667  Mbar,  Poisson's  ratio  of  0.3,  and 
material  density  of  5  gr/cm3 .  Chen  solved  the  problem  using 
the  finite  difference  technique. 

Both  Chen's  solution  and  the  present  solution  are  shown  in 
Figure  4.  A  good  agreement  between  the  two  solutions  was 
obtained . 

As  the  second  example  a  rectangular  plate  with  a  central, 
moving  crack  was  considered  [Ref.  3] .  The  crack  was  assumed  to 
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propagate  at  a  constant  speed  with  uniform  stress  field  of 
68.9  MPa.  The  crack  speed  was  assumed  to  be  33%  of  the 
dilatational  wave  speed  in  a  steel  which  is  5904  m/s.  The 
material  properties  were  elastic  modulus  of  206.7  GPa, 
Poisson's  ratio  of  0.30,  and  material  density  of  8000  kg/m3. 
The  geometry  of  the  problem  is  shown  in  Figure  5.  [Ref.  3] 

A  good  agreement  was  obtained  between  the  analytical  and 
the  present  solutions.  Both  the  analytical  and  present 
numerical  analysis  results  are  shown  in  Figure  6. 

B.  RESULTS  OF  COMPOSITE  PLATES 

A  composite  plate  with  a  central,  moving  crack,  with  a 
similar  configuration  as  the  previous  example,  was 

investigated  for  different  cases.  A  unidirectional  composite 
plate  was  considered. 

For  all  cases  calculations  were  made  for  work  done  on  the 
system,  internal  energy,  kinetic  energy  and  fracture  energy. 
The  objective  of  this  study  was  to  compare  the  fracture 
energies  for  various  conditions.  The  first  case  was  the 
variation  of  the  ratio  of  elastic  moduli  of  the  longitudinal 
and  the  transverse  directions  while  keeping  the  others  fixed. 
The  next  case  changed  cracv  velocities  and,  the  last  case 
varied  material  densities  while  keeping  the  other  parameters 
fixed.  The  geometric  and  material  properties  of  the  problem 
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— 2.4m 
Crack  tip 


Figure  5  Finite  Element  Mesh  for  the  Moving  Crack  Problem 
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are  length  of  3m,  width  of  lm,  elastic  modulus  ratio  in  the 
transverse  direction  of  10  GPa,  Poisson's  ratio  in  the 
transverse  direction  associated  with  loading  in  longitudinal 
direction  of  0.25,  and  uniform  loading  in  the  longitudinal 
direction  of  70  MPa. 

Figure  7  shows  work  and  energy  variations  during  a  crack 
propagation.  For  all  cases  throughout  this  study,  work  was 
higher  than  internal  and  kinetic  energies,  and  internal  energy 
was  higher  than  kinetic  energy. 

In  the  first  case,  the  elastic  modulus  ratios  were 
changed.  Three  different  ratios  between  the  elastic  moduli  in 
the  longitudinal  and  transverse  directions  were  used.  The 
crack  velocity  and  the  material  density  were  kept  constant. 

Figure  8  and  Figure  9  show  work  and  internal  energy 
variation,  respectively.  The  results  reflect  higher  work  and 
higher  internal  energy  for  a  lower  elastic  modulus  ratio.  A 
softer  material  which  has  a  lower  elastic  modulus  ratio  is 
subjected  to  a  larger  displacement  under  the  same  load  while 
all  other  parameters  held  the  same.  This  produces  higher  work 
and  higher  internal  energy  since  they  are  proportional  to 
displacements.  As  will  be  seen  in  all  cases,  the  trends  of 
work  and  internal  energy  are  always  similar  to  each  other. 

Figure  10  shows  the  kinetic  energy  variation.  Kinetic 
energy  is  higher  for  a  lower  elastic  modulus  ratio,  because 
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for  this  material  the  material  velocity  is  higher  due  to  the 
larger  displacement.  Figure  11  shows  the  fracture  energy 
variation.  Fracture  energy  is  also  higher  for  a  lower  elastic 
modulus  ratio. 

The  increments  of  work  and  energies  between  the  initial 
time  and  a  later  time  are  plotted  in  Figure  12.  The  larger 
fracture  energy  for  a  lower  elastic  modulus  ratio  was  caused 
by  a  much  larger  increase  of  work  relative  to  the  increases  of 
internal  and  kinetic  energies  for  this  material. 

In  the  second  case  the  crack  velocity  was  varied.  Three 
different  crack  velocities  were  used.  The  elastic  modulus 
ratio  and  the  material  density  were  kept  constant.  Figure  13 
and  Figure  14  show  work  and  internal  energy  variations, 
respectively.  The  time  needed  for  a  crack  propagation  to  the 
same  crack  size  increment  is  longer  for  a  lower  crack 
velocity.  With  a  longer  time,  the  crack  propagation  influences 
the  rest  of  the  body  more  widely.  The  wider  influence  results 
in  a  larger  displacement.  As  a  result,  it  produces  higher  work 
and  higher  internal  energy.  A  typical  displacement  variation 
at  a  point,  shown  in  Figure  15,  clearly  is  larger  for  a  lower 
crack  velocity. 

Figure  16  shows  the  kinetic  energy  variation.  Kinetic 
energy  is  higher  for  a  higher  crack  velocity.  The  reason  is, 
as  shown  in  Figure  17,  a  higher  crack  velocity  produces  higher 
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speeds  in  the  material.  Figure  18  shows  the  fracture  energy 
variation.  Fracture  energy  is  higher  for  a  lower  crack 
velocity  because  a  lower  crack  velocity  caused  a  larger 
increase  of  work  and  a  less  increase  of  kinetic  energy  than  a 
higher  crack  velocity. 

The  material  density  was  varied  as  the  last  study.  Three 
different  material  densities  were  used.  The  elastic  modulus 
ratio  and  the  crack  velocity  were  kept  constant.  Figure  13  and 
Figure  14  show  the  work  and  internal  energy  variations, 
respectively.  For  the  same  magnitude  of  an  external  force,  a 
smaller  density  causes  a  higher  velocity  and  larger 
displacement  of  the  plate.  Because  the  displacement  is  larger 
for  a  lower  material  density,  work  and  internal  energy  are 
higher.  Figure  22  confirms  it.  Figure  23  also  shows  a  higher 
speed  for  a  lower  material  density. 

Figure  24  shows  the  kinetic  energy  variation.  Kinetic 
energy  is  higher  for  a  higher  material  density.  The  velocity 
is  higher  for  a  lower  density,  but  the  mass  is  larger  for  a 
higher  density.  Because  the  kinetic  energy  is  proportional  to 
the  mass  as  well  as  the  square  of  velocity,  which  case  gives 
a  higher  kinetic  energy  depends  on  which  term  is  dominating. 
In  this  situation  large  mass  yields  a  higher  kinetic  energy 
although  the  material  velocity  is  lower.  Figure  25  shows  the 
fracture  energy  variation.  Fracture  energy  is  higher  for  a 
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lower  material  density  because  of  the  same  reason  as  that  for 
the  case  of  different  crack  velocities.  The  increments  of  work 
and  energies  are  also  plotted  in  Figure  26. 
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Figure  7  Work  and  Energy  Variation  During  a  Crack  Propagation 
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Work  (N.m) 


Figure  8  Work  Variation  for  Different  Elastic  Modulus  Ratios 
and  Fixed  Crack  Velocity  and  Material  Density 
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Figure  11  Fracture  Energy  Variation  for  Different  Elastic 
Modulus  Ratios  and  Fixed  Crack  Velocity  and  Material  Density 
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Figure  13  Work  Variation  for  Different  Crack  Velocities  and 
Fixed  Elastic  Modulus  Ratio  and  Material  Density 
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Figure  14  Internal  Energy  Variation  for  Different  Crack 
Velocities  and  Fixed  Elastic  Modulus  Ratio  and  Material 
Density 
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Figure  16  Kinetic  Energy  Variation  for  Different  Crack 
Velocities  and  Fixed  Elastic  Modulus  Ratio  and  Material 
Density 
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Figure  17  Speed  Values  for  Different  Crack  Velocities  and 
Fixed  Elastic  Modulus  Ratio  and  Material  Density 
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Figure  19  Work  and  Energy  Differences  Between  the  Initial  Time 
and  a  Late  Time  for  Different  Crack  Velocities  and  Fixed 
Elastic  Modulus  Ratio  and  Material  Density 
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Figure  20  Work  Variation  for  Different  Material  Densities  and 
Fixed  Elastic  Modulus  Ratio  and  Crack  Velocity 
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Figure  21  Internal  Energy  Variation  for  Different  Material 
Densities  and  Fixed  Elastic  Modulus  Ratio  and  Crack  Velocity 
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Figure  23  Speed  Values  for  Different  Material  Densities  and 


Fixed  Elastic  Modulus  Ratio  and  Crack  Velocity 
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Figure  24  Kinetic  Energy  Variation  for  Different  Material 
Densities  and  Fixed  Elastic  Modulus  Ratio  and  Crack  Velocity 
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Fracture  Energy  (N.m) 


Figure  25  Fracture  Energy  Variation  for  Different  Material 
Densities  and  Fixed  Elastic  Modulus  Ratio  and  Crack  Velocity 
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Figure  26  Work  and  Energy  Differences  Between  the  Initial  Time 
and  a  Late  Time  for  Different  Material  Densities  and  Fixed 
Elastic  Modulus  Ratio  and  Crack  Velocity 
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IV 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  developed  computer  code  has  been  proven  to  work  well 
for  dynamic  crack  problems .  Very  good  agreements  were  obtained 
for  the  known  solutions.  The  parametric  study  of  dynamic  crack 
propagation  in  composite  plates  provided  the  following 
conclusions . 

A  softer  composite  plate  which  has  a  lower  elastic  modulus 
ratio  showed  a  higher  fracture  energy  if  the  crack  velocity, 
the  material  density  and  other  characteristics  of  the  problem 
were  held  constant.  A  composite  plate  with  a  lower  crack 
velocity  and  a  composite  plate  with  a  lower  density  showed  a 
higher  fracture  energy  provided  all  other  parameters  were  held 
constant,  respectively. 

In  this  study  a  discontinuous  crack  release  technique  was 
used  to  model  the  crack  propagation.  This  causes  sudden  shock 
waves  which  in  turn  result  in  oscillations  in  solutions.  A 
technique  that  releases  the  propagating  crack  nodes 
continuously,  for  example  as  shown  in  Ref.  12,  is  recommended 
to  achieve  more  stable  solutions. 
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APPENDIX 


The  general  stiffness  matrix  using  a  bilinear  rectangular 
element  and  the  general  case  of  material  property  [D]  matrix, 
that  is. 


is  given  as  below: 


[D]  = 


*11 

*12 

*13 

*21 

*22 

*23 

*31 

*32 

*33 

*11 

*12 

*13 

*14 

*15 

*16 

*17 

*18 

*21 

*^2 

*23 

*24 

*25 

*26 

*27 

*28 

*31 

*32 

*33 

*34 

*35 

*36 

*37 

*38 

*41 

*42 

*43 

*44 

*45 

*46 

*47 

*48 

*51 

*52 

*53 

*54 

*55 

*56 

*57 

*58 

*61 

*62 

*63 

*64 

*65 

*66 

*67 

*68 

*71 

*72 

*73 

*74 

*75 

*76 

*77 

*78 

*81 

*82 

*83 

*84 

*85 

*86 

*87 

*88. 

where 


b2c2  +  M-D33b*C 

K12  =  lfD12bc^*  {D12+DZ3)  b2c2  +  lfD32bic 
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J<ri3=-^DIlJbc3+4{-D31+D13)i2C2  +  |D33Jb3C 
K1A=-±$-Dx3bc*+4  iD12-D,2)  b2C2  +  ^D32b*C 
K1S  = - 1  D^bc 3  +4  ( -D3 2 -d13 )  b 2 c 2 - -| D^b 3 c 
K16  =  -^D13bc^A(-D12-DZ2)b2C2-^DZ2b2C 
7  =  |  i  ^  3  +  4  ( D3  ,  -  Dx  3 )  2) 2  C  2  -  D3  3  2) 3  c 

Kls  =  ^D13bc2+4  (~Dl2+DZ3)  b2c2-^DZ2b2c 
%2i= ~y~D31bc3 +4  (021+D33)  b2 cz + ■—■  D23b2 c 
^^D^bc2**  (DZ2+D23)  b2C2  +  ^D22b2c 
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^-^D^bc'+i  (~D21+D22)  b2c2  +  ^D22b2C 
K2i  =  -^fD22bc2^i  {D22~D22)b2c2  +  ^D22b2c 

K25  =  -^D31bc2+U-D21-DZ3)  b2c2-^D22b2C 


K26  =  -^D22bc2+4(~D22-D22)b2C2-^D22b2C 


K,  ,  =  |  D2  x  b  C 3  +  4  { D2  2  -  D2  3 )  b 2  c  2  - D2  3  b 3  c 


K2S  =  -^D22bc2+4  (~D22+D22)  b2C2-^j-D22b2C 


K21  =  -^D21bc2+4  {D21-D12)  b2C2  +  ^D22b2C 


K22  =  ^fD12bc2+U-Dl2+D22)b2c2  +  ^D22b2C 
^  —  D^bc2^  (-D21-D12)  b2c2^D23b2C 
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*34=^13^3+4(  -^12-^33)  b2C2  +  ^-D32b3C 

K25  =  —  D11bc2  +4  (~D3X+DX3)  b2 C2  -  -j- D33b3 C 
K3€=^Dx3bc3^  ^x2-D33)  b2  c2  -  D32b3  c 

K21  =  -^Dxxbc3^  (D3X+DX3)  b2C2-^D33b3C 
K3^-^Dx3bc3^(Dx2*D33)b2C2-^D32b3C 
K*x  =  ~D3xbc3+l{D2X-D33)b2C2  +  ^D23b3c 
K*2~~ D33bc3  +4  ( -D32+D23)b2c2+^D22b3c 


K<2  =  Jj-D3xbc3+4(-D2X-D33)b2C2  +  ^-D23b3C 


63 


K*4  =  JyD23bci+4  (-£32-D23)  b2C2  +  ±^D22b2c 

*4 5  =  - f D2 V bl C' 3  +  ' 4  (  -■ °2 1 +1 °3 3  )  b“ 2‘ C' 2  -  D2 3i b  3 C 
Ki6  =  ^D32bc2^A  (D22~D23)  b2c2~^-D22b2C 
K„  =  -^D3xbc2+±[D21+D33)b2c2-^D23b2C 
K^-^D33bc2*A  (D32+D23)  b2C2-^D22b2C 
K5^-^Dxlbc2*A(-D3X-Dx3)b2C2-^D33b2C 
Ks2  =  ~^Dx3bc2+A  (-DX2~D33)  b2C2-^D32b2C 
K53=^Dxxbc^A  (D3X-DX3)  b2C2-^-D33b2C 
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Ks4  =  ^Di3bc2+U-D12+D23)b2C2-l£-D32b2c 


*55  =  D^bc*  +4  { 03l +D12 )  b: 2 C: 2  +  D33 b 3 C 

JC56  =  -y-Di3iDC3+4  (D^D^)  b2C2*~DZ2b3C 

*5,  =-^j-Dlrbc3+4  ( -d31+d13) b2c2  +  ±$-D33b2c 
Kss  =  -^fD13bc^A(D12-D33)b2C2  +  ^D32b2C 

*6i  =  -|^3i^3+4(-D21-i?33)Jb2C2-|D23Jb3c 

^2  =  -|^33^3+4(--D32--D23)^2^2-|I)22^3C 
^  3  =  |  d3  ,  i3  C  3  +  4  (  d2  ,  -  z:3  3  )  i  2  C  2 d2  3  i3  3  C 
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%65  =  D3lbc2  +4  (D21+D2^)  b2  c2  +  -^-  D23b2  c 
^^J-V^bc2--*  (D32+D23) b2c2+^D22b2c 

K61=-~-D31bc2+4  (-D21+D33) b2c2+^j-D23b2c 

K^8  =  -^-D33bc2*4  (D32~D23)  b2C2  +  ^D22b2C 

K, t  =  £ D13bc 2  +4  (  -D31  +D13 )  b' 2C2  - D33b 2 C 

K,2  =  ^D12bc2+4(D12-D33)b2C2-lfD32b2C 

JC73  =  --|D11i3C3+4  (D31+D13)b2C2~^D33b2C 
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J^4  =  - 1  D13i)C 3 +4  (D12 +D33 )  i 2 C 2  - 1  D32,b 3  c 


%75  =  ~~:j~D11bc2  +4  (D31-D13 )  b2c2  + -^ D33b2 c 


*,6 =' -• ~ D^c 3 +4  < -^12 +^33 )  b 2 c2 * | D32b 3 c 


K71  =  ^j-D11bc2^(-D31-D13)b2C2^D33b2c 


^S=~D13bc2^(-D12-D^)b2C2*^-D22b^c 


K81  =  ^D31bc^A(-D21+D33)b2c2-^-D23b2c 


KS2  =  ^D33bc^4  (D32-D23)  b2c2-^-D22b2C 


K93  =  --^D3lbc2+4(D21+D33)b2C2-£D23b2C 
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8 


K.A  =  -4D„bc3+  4  (ZX,+Z?,,)2?2C2-- |s,,h3C 


v84 


33 


32  ^23 > 


22 


^85“  — 2~^31^>C  3  ^21  ^33^  b2C2  +  —  D22b2C 


K86  =  -lf-D33bc>+U-D22+D23)b2C2  +  £D22b*c 


K81~ D3lbc3 +4  (  D2X  D33)b2C2  +  ^D23b3C 


K88  =  lj-D33bc3+4(-D32-D23)b2C2  +  ^-D22b3c 


In  this  study  the  [D]  matrix  for  a  unidirectional  composite 
material  for  the  plane  stress  conditions  is; 


[D] 


D. 


21 


012 

022 

0 
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and  the  nonzero  elements  of  the  [D]  matrix  are; 


Dll= 


12  21 


^22  ”  ' 


i-*  12*21 


n  _n  _  U12^2  _  U21^1 
-^12 


1  U12U21  ^  U12U21 


^33  ^12 


where  subscripts  1  and  2  on  the  right  hand  side  of  the 
equations  are  the  direction  of  the  fibers  and  the  transverse 
direction,  respectively.  Substituting  these  values  into  the 
general  stiffness  matrix  gives  the  stiffness  matrix  for  a 
unidirectional  composite. 
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